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ABSTRACT 



J3 • 

We demonstrate an approach to solving the coagulation equation that involves 

using a finite number of moments of the particle size distribution. This approach 

^j ! is particularly useful when only general properties of the distribution, and their 

time evolution, are needed. The numerical solution to the integro-differential 

Smoluchowski coagulation equation at every time step, for every particle size, 

and at every spatial location is computationally expensive, and serves as the 

Q\ \ primary bottleneck in running evolutionary models over long periods of time. The 

00 

advantage of using the moments method comes in the computational time savings 

gained from only tracking the time rate of change of the moments, as opposed to 

tracking the entire mass histogram which can contain hundreds or thousands of 

bins depending on the desired accuracy. The collision kernels of the coagulation 
O 



equation contain all the necessary information about particle relative velocities, 
cross-sections, and sticking coefficients. We show how arbitrary collision kernels 
may be treated. We discuss particle relative velocities in both turbulent and non- 
turbulent regimes. We present examples of this approach that utilize different 
collision kernels and find good agreement between the moment solutions and the 
moments as calculated from direct integration of the coagulation equation. As 
practical applications, we demonstrate how the moments method can be used to 
track the evolving opacity, and also indicate how one may incorporate porous 
particles. 

Subject headings: accretion, accretion disks — methods: numerical -- radiative 
transfer — solar system: formation 
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Introduction 



The conditions under which planetary formation is initiated in the protoplanetary disk 
remain poorly understood despite several decades of astrophysical research. Furthermore, 
the continued discovery of extrasolar planets (over two hundred currently) further em- 
phasizes the need to understand the formation of not only our own solar system, but of 
these numerous strange and diverse systems as well. At the forefront still remains the 
basic question of how growth occurs from dust to planet. The key properties of proto- 
planetary nebulae remain controversial, yet grain growth from dust to larger agglomerates 
(~ cm-m), and from planetesimals (> km) to planets must have occurred in some man- 
ner. While numerous studies of iV-body dynamics have addressed the growth of terres- 
trial planets from asteroid-size planetesimals which are l argely decoupled from the gas (e.g., 
Leinhardt fc Richardsonll2005l ; iBromley &: Kenyonll2006l ). it is the transition from agglomer- 
ates to planetesimals which continues to provide the major stumbling block in planetary ori- 



gins (IWeidenschi lling 



1997 



2002 



2004; IWeidenschilling fc Cuzzilll993l; Istepinski & Valageas 



6: 



19971 ; iDullemond &; Dominiki 120051 ; ICuzzi fc Weidenschilling 



200 



Dominik et al. 



2007m . 



Grain growth begins with sticking of sub-micron sized grains, which are dynamically 
coupled to the nebula gas and collide at low relative velocities which are size-dependent 
(jWeidenschilling fc Cuzzil 119931 ). These relative velocities can be caused by a variety of 
mechanisms, such as Brown ian motion for the smallest grains, differences in coupling to 



eddies in a turbulent nebula ( Volk et al 



198 



0; 



Weidenschilling 



1984 



Markiewicz et al. 1991 



Weidenschilling fc Cuzzil Il993l ; ICuzzi &: Hogan 



2003 



Ormel fc Cuzzi 



radial, and azimutha l drift driven by global gas pressure gradients (INakagawa et al.lll986 



2007) , and/or vertical 



Weidenschillingi 119971 ). For the smallest sized grains and aggregates, sticking is caused by 
weak van der Waal s interaction forces (although stronger forces may also act; see, e.g., 
Dominik et al.l 120071 ). This sticking forms larger and larger aggregates until particles grow 
large enough to decouple from t he gas and settle to th e midplane where they may, barring 
strong turbulence ( althoug h see iJohansen et al.l 120071. for an exception) , grow into larger 



objects, even planetes imals (jSafronov 
Weidenschillingi 1 19971 ) . 



1969 



Weidenschilling fc Cuzzilll993l ; ICuzzi et al.lll993 



There is considerable observational evidence supporting dust growth in protoplanetary 
disks, at least from sub-micron to millimeter scales, based on observations of the dust con- 
tinuum emission originating from protoplanetary disks. The wavelength dependence of dust 
emission from IR to radio wavelengths suggests, in many cases, that typical grain sizes 
are in the millimeter rang e , which is orders of magnitude larger than inters tellar grains 
(JBeckwith fc Sargentl 1 199 It IDullemond fc Dominiki 120051 ; iDominik et al.l 120071 ) . Measure- 
ments of the contrast of the Orion trapezium silicate emission feature at 10 /im relative to the 
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local continuum implies an increase in size from interstellar dust (sub-m i cron) to micron sizes 



in th e disk photospheres (Ivan Boekel et al.ll2003l ; iPrzygodda et al.ll2003l ; iKessler-Silacci et al. 
2006h . 



Theoretical approaches to the study of dust coagulation normally involve solving the 
cumbersome collisional coagulation equation (see Eq. 1) at each spatial location (R, 0, Z) 
in the disk. The collisional kernel K(m, m!) for particles of mass m and m' usually in- 
volves some sort of sticking efficiency, and particle cross-sectional area, in addition to the 
particle-particle relative velocities. Although the systema tic particle-particle velocities that 
arise from gas pressure gradients (INakagawa et al.l Il986l ) are somewhat easy to incorpo- 
rate into the collisional kernel, relative velocities in turbulence are mor e complicated. I f 
one assumes that nebula turbulence follows a Kolmogorov inertial range (IVolk et al.lll980l ) 
in which energy cascades from the largest, slowly rotating eddies, down to some minimum 
lengthscale where the intrinsic molecular viscosity can dissipate the macroscopic gas motions 
(and turbulence ceases), one can obtain c losed-form expressions for the particle-particle, and 
particle-g as turbulent relative velocities (ICuzzi fc Hoganl 120031 ). even for particles of differ- 
ent sizes (jOrmel fc Cuzzil 120071 ). Numerical models of particle growth in the inner regions of 
protopla netary disks indicate that incl usion of these turbulence-induced relative velocities is 
essential (JDullemond fc Dominikil2005l ). as it provides a source of energetic collisions that, in 
the presence of fragmentation and destruction, leads to the ongoing existence of small grains 
even after several million years. 

A full scale solution to the problem of dust coagulation for a size spectrum at every 
vertical and radial location in the protoplanetary disk, as a function of time, is th us pro- 
hibitive even with today's computational capabilities (jWeidenschillingl 120021 . 12004 ) . So it 
is advantageous to find alternative methods of extracting information from the coagulation 
process, such as time scales of growth, the largest particle size, the size that dominates the 
area (radiative transfer), and the size that dominates the mass (migration and redistribution) 
without having to track the behavior of the entire mass distribution. 

In this paper, we describe a method that utilizes the moments of the coagulation equa- 
tion to accomplish this task, in particular when the functional form of the mass distribution 
can be assumed. In § 2, we derive the moment equations from the coagulation equation, and 
compare direct integration of the coagulation equation with solutions using the moments 
approach in the case of a simple turbulent coagulation kernel. We also present two different 
approaches to obtaining solutions using a realistic kernel, when the form of the size distri- 
bution is assumed to be a powerlaw. In § 3, we evaluate the different approaches developed 
in § 2 for realistic kernels by numerical integration of the moments equations and compare 
these moments to those obtained by direct integration of the coagulation equation. We also 
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compare the moments solution to an alternate analytical approach by iGaraudl (120071 ) . In § 
4, as a demonstration of the practical application of the moments method, we calculate the 
wavelength-independent opacity, and discuss generalization to wavelength-dependent Planck 
or Rosseland mean opacities. In § 5, we indicate how porosity may be included. In § 6, we 
present our conclusions. 



2. Solution Techniques Using the Moments Method 

The moments method approach to modeling particle growth allows an attractive alter- 
native to solving the full coagulation equation when it is only necessary to know a particular 
property, or propert ies of the evolving particle size distribution. The coagulation equation 
(jSmoluchowskil Il916l ) in its integro-differential form is given by 



df{m,t) 



dt 



K(m — m' , m')f(m — m', t)f(m', t)dm! 



in 

Kim, m')f(m, t)f(m', t)dm' + — 

o dt 



d£ 
dt 



where f(m,t) is the particle number density per unit mass at mass m, and K(m,m') is the 
collision kernel connecting the properties of interacting masses m and mf, which typically 
takes into account mutual particle cross section a(m,m'), relative velocity AV(m, mf), and 
a particle sticking efficiency S{m,m'). The last two terms on the RHS of Eq. (1) represent 
sources and sinks such as particle erosion, fragmentation and subsequent redistribution of the 
fragmented population, and gravitational growth. Although we will not specifically address 
these mechanisms, we will indicate how they may be treated, and save implementation for a 
forthcoming paper. 

The motivation behind using the moments method is to avoid the computational cost 
inherent in Eq. (1) which entails solving the convolution integral (the first term on the RHS 
of Eq. 1) at every location, and at every time step. Depending on the functional form of the 
kernel, the second integral on the RHS may also need to be integrated repeatedly. Given that 
the typical mass spectrum may involve 10 2 — 10 3 bins to acquire the desired accuracy, the 
computational burden required becomes a detriment to any study involving a wide range 
of parameter space. The so-called brute force solution of the coagulation equation thus 
becomes the primary bottleneck in running 2D evolutionar y models over extended periods 
of time (jWeidenschillingl 120041 ; iDullemond fc Dominiki 120051 ) . 



We define the p-th moment M„ of the distribution, where p need not be an integer, as 



M n 
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m p f(m,t)dm, 



(2) 



where the units of f(m,t) are such that Mq = J f(m,t) dm represents the total number 
density of particles, and M± = J mf(m, t) dm = p is the total volume mass density of solids. 
The essence of the moments method is as follows. We multiply both sides of the coagulation 
equation (Eq. 1) by m k , where k is an integer, and then integrate both sides over mass m: 



dM k 
dt 



kdf, 1 

m — dm = - 
dt 2 



m h dm / K{m — m',m')f(m — m',t)f(m',t)dm' 



OO fOO 

k . 



m f(m,t)dm / K(m,m')f(m',t)dm'. (3) 
i) Jo 



We introduce a step function H(m — m'), such that H = for m — m! < and H = 1 
otherwise, to extend the limits of the integral over m! from (0, m) to (0, oo). The convolution 
integral (the first integral on the RHS of Eq. 3) then becomes 



-1 /*oo /*oo 

- / m h dm / H(m — m')K(m — m',m')f(m — m',t)f(m',t) dm' 
2 Jo Jo 



, f(m',t)dm' \ (u + m) h K(u,m')f(u,t) du, 
2 Jo Jo 



(4) 



where on the RHS of Eq. (4) we have switched the order of integration and made the 
substitution u = m — m! . The RHS side of Eq. (4) may then be combined with the last 
double integral in Eq. (3) to yield the set of ord inary differential equations (ODEs) for the 
integer moments (IMarov fc Kolesnichenkoll200ll ): 



dM k 
dt 



OO /"OO 



JO 



1 



m + m') k - m k 



K(m, m')f(m, t)f(m', t)dmdm', 



(5) 



where we have substituted m for u with no loss of generality. Depending on the mass 
dependence of the kernel (as illustrated below) the right hand side can readily be expressed 
as products of moments of order k or less, leading to a closed system of equations which may 
be solved using standard techniques. In this definition of the coagulation equation in which 
it is assumed there are no sources and sinks, the first moment Mi = p is constant in time, 
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that is dMx/dt = C0. 



Exact solutions to the coagulation e q uation have been obtained for some specific choices 



of the kernel (e.g., ISmoluchowskil Il916l ; iTrubnikovl Il97ll ). the most simple being that of 



constant kernel K(m,m') = (3q. Although the exact solution for f(m,t) cannot be obtained 
from the moments equations, the time rate of change of the zeroth moment (k = 0), can 
easily be obtained from Eg. (5) which red uces to dM /dt = — (1/2)/3 Mq and has the trivial 
solution (see, e.g.. ISilk fc Takahashilll979l ) 



M (t) 



M o (0) 



(l/2)/3oM (0)f 



(6) 



independent of the initial choice of distribution f(m,0). Likewise, the ODE for M2, which 
can be associated with the density- weighted mean particle size ((m) = M 2 /Mi), yields the 
simple solution M 2 (t) = M 2 (0) + j3 Q p 2 t. Despite not knowing the exact solution, we are able 
to understand the behavior of general properties of the mass distribution with time through 
the moments equation without tracking the behavior of the full mass spectrum. Thus, if it 
is only desired to know, for example, the time evolution of the particle representing most of 
the area (first moment) or most of the mass (second moment) of the distribution, then the 
advantage of the moments method becomes clear. A small number of moments is all that is 
necessary to determine the behavior of the system. In particular, we will be interested in the 
size of the largest particle rriL{t) in the entire mass distribution as a function of time (which 
we show below can be computed from the integer moments ), because these are usually th e 
most rapidly drifting and most violently colliding particles ( jCuzzi fc Weidenschillingll2006l ). 



2.1. Example: Saffman and Turner Turbulent Coagulation Kernel 



We can illustrate the moments method appro ach using a very simple turbulent coagu- 
lation kernel where K(m,m') = 7o(m 1 ' 3 -l-m' 1 ' 3 ) 3 ( jSaffman fc Turnerlll956l ). and no sources 
or sinks. Physically, this represents the product of (r + r') 2 for area, and (r + r') for relative 
velocity of two particles in a laminar shear flowfield. If a sticking coefficient were desired then 
we would specifically have K(m,m') = 7o(m 1//3 + m n ^ 3 ) 3 S(m,m') where < S(m,m') < 1. 
Here, 70 is a constant that depends on the Reynolds number of the gas flow. Inserting this 
kernel into Eq. (5), we find the set of equations 



1 We note that under realistic protoplanetary nebula conditions, p will not be constant due to, e.g., size- 
dependent advection terms in the equations. Such effects can be treated separately from the "coagulation" 
step. See § 2.2.2. 
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dM 
dt 



- l0 (M M 1 + 3M 1/3 M 2/3 ), 



dM x 
dt 



dM 2 
~df 



2 7 o(M 2 M 1 + 3M4/3M 5/3 ). 



(7) 



We note that the physical derivation of kernels in terms of particle radius r leads to fractional 
moments in terms of m. These fractional moments look complicated, but can be solved for 
by sim ple interpolation using Lagrange polynomials in terms of the mo re familiar integer mo- 
ments ( )Loginovill979l ; iPress et al.lll992t iMarov fc Kolesnichenkoll200ll ). Thus, any fractional 
moment M p may be expressed compactly in the normalized form M' At) = M p (t)/M p (0): 



k+n 

kw = n [Mi®] 

j=k 



L"(p) 



where n is the number of integer moments, k = 0, 1, ..., n, and the exponent LP is defined as 



TV. J -n+l 



P-3 



(9) 



with EL+iG ) = Pip ~ 1) — (P — n )i an d the n are the binomial coefficients n\/j\(n — j)\. 
The important thing to note here is that the order of the moments must remain less than or 
equal to the largest moment in order for the system to remain closed. In general, this will 
be true for realistic collision kernels (see § 2.2). 

To test the accuracy of the moments method, we integrated equations (7) using a fourth 
order Runge-Kutta scheme, and compared the results to a brute force integration of the co- 
agulation equation (Eq. 1). For the latter, we integrated the distribution function f(m,t) 
at each timestep to determine the moments of the distribution Mo, Mi, and M 2 . We chose 
f(m, 0) = CQin~ q as our initial distribution for simplicity. Since the form of the mass distri- 
bution is only assumed at t = 0, we consider this to be an example of an implicit assumption 
(see § 2.2.2). The results of this calculation are presented in Figure 1. We have used two 
different resolutions for the brute force calculation in order to demonstrate how the higher 
resolution (and much more computationally expensive) case (solid symbols) approaches the 
moments approach solution. Notice that the first moment Mi = p remains constant as ex- 
pected. Given that general numerical errors can arise from the finite mass grid and coarse 
timesteps (At = 10 years) used for the calculation, and that systematic errors may be in- 
troduced by the Lagrangian interpolation, the fit is quite good. It is important to point 
out here that while the coagulation calculation required as many as 20 — 30 hours of CPU 



time to complete on a 2 GHz machine, the moments calculation of Eq. (7) is essentially 
instantaneous. 



2.2. Realistic Coagulation Kernels 

The Saff man- Turner turbulent coagulation kernel we used as an example in § 2.1 is 
simple to utilize and is expressible explicitly in powers of the mass m. In practice, however, 
the realistic coagulation kernels that we will be using will be more complicated than this 
simple example. A realistic coagulation kernel will be, at the very least, a product of a 
mutual cross section 7r(r + r') 2 and a relative velocity 



AV(m, m') 



(U - U') 2 + (V - V') 2 + (W - W') 2 + v 2 



m,m 



(10) 



where U(m), V(m), and W(m) are the x,y,z systematic (pressure gradient driven) particle 
velocities, and v? is the turbule nt velocity contribution which in general is n ot separable into 
distinct functions of m and m! (ICuzzi fc Hoganll2003l ; lOrmel fc Cuzzill2007l ). The systematic 
velocities, which are derived for a discretized particle size distribution in the appendix, are 
complicated functions of the particle size through the stopping time t s , 



mAVr 



i>u 



;n) 



D 



where AV pg = |V — v| is the relative velocity between the particle and the gas, and Fd is 
the drag force on the particle of size r which dep ends on the size o f the particle relative to 
the mean free path A of the gas molecule (e.g., see lCuzzi et al.lll993l ). Depending on whether 
r > A or r < A determines whether the stopping time itself depends on the instantaneous 
relative velocity AV pg between the particle and the gas (Stokes regime) or does not (Epstein 
regime). In the latter case, calculation of the systematic velocities (U, V, W) and gas velocities 
(u, v, w) for a particle size distribution is straightforward. In the former case, iterations are 
required to correctly calculate AV pg (see appendix). 

The turbulent velocities are less straightforward to implement than their systematic 
counterparts because of the different coupling that exists between particles and eddies of 
different sizes (and thus stopping and decay times, respectively). It has not been until very 
recently that closed form expre ssions for the turbule nce-induced velocities (for a particle 
size distribution) were derived (jOrmel fc Cuzzil 120071 ) . These expressions can be written 
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separately for the so-called "class I" and "class II" eddiesj as: 



1 



St 



St 



st? 



+ 



st? 



Sti + St,- V St*, + Sti 1 + Sti St*, + St, 1 + St 



(12) 



AV 2 T 



where v\ 



2(St£ - Re~ 1/2 ) 



Sti 



St 



Otj + i^tj.- 



St, + Re" 1/2 



+ 



st? 



St 



Dt^' ~r O L^ y 



St, + Re- 1/2 



(13) 



AV} 2 + AV}j, Re i s the Reynolds number , v g 



a 



1/2 



c 9 is the turbulent ga s velocity 



with c 9 the sound speed (e.g-. lNakagawa et al.lll986l ; ICuzzi fc Weidenschillingll2006l ). and the 
particle Stokes numbers Stj = t^/tz, where th is the turnover time of the largest eddy, 
typically taken to be ft' 1 . The boundary between class I and class II is defined by the 
"combined" Stokes number St*, =_max(St*, St*) and the values of the St*, are obtained from 



the equation (jOrmel fc Cuzzil 120071 ) : 



•-.! 



:Vk(yl 



St, 



1+tf 



1 + Sti 



+ 



1 AV 2 
1 pg 

St*, u 2 



(14) 



where y| = St*,/Stfc. In our calculations, we will make use of these expressions when com- 
paring cases with turbulence-induced velocities. 

Given that the moments method requires that we express the kernel in terms of fractional 
or integer moments, the problem becomes expressing the relative velocity (Eq. 10) in a form 
that satisfies this requirement. As an example, the systematic velocities U, V, and W can 
each be fit easily by a finite series in fractional powers pi of m as U(m) = Yli CLi^n Pl where 
the coefficients o, can be found by fitting N points to the function U(mi) = Ui (I an integer). 
This leads to the system of equations 



N 



« = £ 



aim 1 



(15) 



i=l 



2 The concept of eddy classes were introduced by IVolk et al.l (jl98df ) to distinguish between slowly and 
rapidly varying eddies. Class I eddies are defined as those in which the eddy fluctuates slowly enough that a 
particle's stopping time t s is much shorter than the eddy decay time and the time it takes to cross the eddy; 
thus, they will align themselves with the gas motions of the eddy prior to its decay. Class II eddies then are 
defined as ones in which the eddy decay time is fast compared to the particles t s , and thus only provide a 
small perturbation to the particle motion. 
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The coefficients a, then follow from matrix inversion. Similarly, one can find expressions 
V(m) = Y2i biin Pi , and W(m) = ^ Ciirt Pi so that we may construct the laminar expression 
for (AV) 2 in terms of the variables m and m! by multiplying out the individual terms, e.g., 
U(m)U(m') = Y^i A(ii(ijV(i Vi m! Vi and so on. We find then that we can express the square of 
Eq. (10) as 

(AV(m, m')) 2 = J2 A ij [m Pl+P3 - 2m p >m' p > + m' p '+^] , (16) 

hj 

where A^ = a^aj + b^bj + CiCj, and Pi + pj < 2 for closure. This determines the matrix of 
coefficients A^ which operate on the finite power series in m and m'. Although we cannot use 
the same approach for solving for AV in Eq. (10) because of the radical, we were motivated to 
express AV in a similar form, i.e., as AV(m, m!) = ^ • A' i Am Pi+Pj — 2m Pt m' Pj +m' Pi+Pj ], and 
Pi+Pj < 1- The problem reduces to solving for the coefficients A'- directly, by similar matrix 
inversion techniques. Unfortunately, the number of points needed to obtain an accurate 
representation of the two-variable function AV(m, m') this way exceeds the limitations of 
the inversion. In practice we found that we could solve a system for AV with at most 
5 — 6 points (25 — 36 matrix elements), whereas Eq. (16) carries much less restriction in 
the number of points needed because we could construct (AV) 2 by the product of accurate 
representations of single variable functions. Furthermore, due to the coupling evident in 
the turbulence induced velocities (Equations 12 and 13), the turbulent velocities are not 
separable functions of the masses m and m' , further complicating the analysis. 

Powerlaw assumption: Fortunately, it turns out that, although the direct calculation 
approach to AV described above would be mathematically appealing in the sense that the 
moment equations would remain implicit (i.e., the form of the mass distribution would only 
be assumed at t = 0), it is not necessary. In the next sections we describe two alternative 
approaches to dealing with realistic coagulation kernels that employ the moments method 
under the assumption that the form of the mass distribution f(m, t) is a powerlaw. Such an 



assumption is not entirely unfounded. A number of detailed models by lWeidenschillingi (11997 



20001 . 12004 ) have shown that powerlaw size distributions result, which have nearly constant 
mass per decade radius to an upper limit mx,(t) which grows with time until the frustration 
l imit of around a meter is reac hed, and (under turbulent con ditions, at least) growth stalls 



( jCuzzi fc Weidenschillingl 120061 ) . Similar trends are found by iDullemond fc Dominik! (120051 ) 
in which different assumptions about collisional ejecta are made. These authors do find 
minor fluctuations in the distribution, but if one is primarily interested in general properties 
of the distribution, such as how the largest particle size changes with time, and not the 
fine structure of the mass distribution itself, the approach is quite advantageous. There are 
reasons to believe that a powerlaw is a natural end-state, especially those with equal mass 
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per decade, because they have self-preserving properties for the collision kernels of interest 



( ICuzzi fc Weidenschillind 120061 ) 



The significance of the upper mass cutoff m^ depends on the assumed slope of the 
powerlaw. In all real distributions, there will be a rapidly decreasing abundance of particles 
for masses exceeding some threshold, even though the abundance may not drop immediately 
to zero as in our assumed model. The rare, extremely large particles might be of interest for 
some applications, but our focus will be the particle size carrying most of the area, or most 
of the mass (the first or second moments). For powerlaw distributions f(m) oc m~ q with 
q < 2, rriL is itself the mass of the particle carrying most of the mass and is independent 
of the selection of a lower particle size cutoff. For q = 2, the distribution contains equal 
mass per decade, and for q > 2, most of the mass is in the small particles and the value 
of rriL depends on the lower particle size chosen. Most realistic distributions, and those 
most commonly treated in the literature, have q < 2; here we assume q = 11/6, a widely 
used fragmentation powerlaw. In this case, a strict cutoff at m L represents a well-defined 
distribution with easily-understood moments where most of the mass is at tul and the area 
is nearly equally distributed per decade with a mass dependence vrC 1 ^ . 

Although we do not include either imperfect sticking or fragmentation at this stage in the 
model, we believe that the moment equations as expressed in Eq. (5) will remain valid up to 
a "fragmentation barrier", which may be defined as that size for which the typical disruption 
energy of a particle is on the same order as the energy of identical colliding particles. The 
fragmentation barrier will also depend on one's choice of nebula parameters. T h is tre atment 



(up to the fragmentation size) is consistent with recent work by IBrauer et al.l (120081 ) (their 
Fig. 13) which shows a constant powerlaw mass distribution up to a cutoff size which then 
falls off abruptly. This "knee" in the distribution represents that efficient fragmentation size. 

Once mi reaches the fragmentation barrier, growth beyond this stage would need to 
be treated in a different manner, for exam ple, by simple sweepup of small, less disruptive 



particles by large ones (ICuzzi et al.lll993l ). Creation and disruption of these particles can 
be handled as part of the source and sink terms in which for example, disrupted particles 
are assumed to be fragmented back into a powerlaw distribution (which is suggested by 
experimental evidence and widely assumed by other modelers) as opposed to monomers. 
The model as presented within this paper, however, should be useful for the early stages 
of protoplanetary nebula particle growth relevant to spectral energy distributions and MRI 
suppression. We leave the incorporation of growth stages beyond the efficient fragmentation 
stage for a later paper. 
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2.2.1. Approach 1: Explicit Assumption 

Motivated by our discussion of the previous section, if we assume that the form of 
the particle mass distribution function remains a powerlaw at all times, we may express 
f(m,t) = c(t)m~ q , such that m L {t) is the growing upper limit of the distribution, q is the 
slope which is assumed to be constant (although see below), and c(t) is a normalization 
coefficient. Taking the lower limit of the mass distribution to be mo, then moments as 
expressed in Eq. (2) are explicitly given for q — p ^ 1 by 

M p (t) = - C(t) {m L {t) l+p - q ~ ml +p - q ) . (17) 

1 + p — q 

We then take the time derivative of Eq. (17) for the zeroth and second moments, and 
substitute these expressions on the LHS of the corresponding moment equations in Eq. (5) 
to get 

where Eq. (18) is valid for q ^ 1, and Tq and T 2 are the integrals on the RHS of Eq. (5) for 
k = and k = 2, respectively. That is 

rm L (t) rm L (t) 

Tq^ttil) = / / K(m,m')m~ q m'~ q dm dm', (20) 

J mo J mo 

rm L (t) rm L {t) 

T 2 (m L )= I / K{m,m')m l - q m' l - q dmdm' ) (21) 

where the kernel is given by K(m,m') = o~(m,m')AV(m,m')S(m,m'), with a(m,m') = 
Ko{m 1 ^ + m' 1 / 3 ) 2 , Kq = 7r(3/47rp s ) 2 / 3 , and p s is the particle material density, assumed 
constant. Note that for the special case of q = 1, Eq. (18) has a slightly different form which 
depends on In {m^jm^). After some simple algebra, Eq.'s (18) and (19) can be written as 



drriL 

— ; — = c 
dt 



(3 - q){m l L q - m\ q )T 2 + f (1 - q)(m 2 L q - m\ q )T 
(3 — q)(m L ~ q — m^ qs )m^ q — (1 — q)(m L q — m ~ q )m~ L q 



(22) 
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dc 9 

dt 



(3-g)(l-g)(m L 9 r 2 + i 



m 



2-q, 



(3 — q){m L q — m q )m L q — (1 — q)(m L q — m q )m 



(23) 



which we integrate using a fourth order Runge-Kutta scheme. Equation (17) then gives the 
moments Mq and M 2 as a function of time, which can be directly compared with direct 
integration of the same conditions using the coagulation equation (Eq. 1). 

The advantage of this approach (in which a powerlaw is assumed at all times) is its 
transparency; that is, the variables being sought (mi and c) are solved for directly. Fur- 
thermore, the change in the coagulation kernel as the particle size distribution changes is 
included because the kernel is updated and explicitly integrated into Tq and r 2 with every 
time step. This will prove advantageous when additional effects such as sticking are included 
in the kernel. In addition, source and sink terms need not parameterized in terms of integer 
moments and can be implemented directly. An unfortunate disadvantage of this approach 
is that because the kernel must be integrated (in fact several times) over both m and m! 
every time step to get To and r 2 , the CPU time involved is significantly longer than a fully 
implicit case (e.g., § 2.1; also see § 2.2.2); however, it remains a much faster approach (or- 
ders of magnitude) than solving Eq. (1) directly since the cumbersome convolution has been 
eliminated. 



2.2.2. Approach 2: Semi-implicit Assumption 

Here, we present an alternative moment-based approach to solving the realistic coag- 
ulation kernel in which the integer moments appear directly. Unlike the previous case of § 
2.2.1, where the double integral of Eq. (5) is used to get the functions of m L , T and T 2 , we 
only integrate over one mass variable (i.e., only integrate one of the integrals), defining the 
functions 

C (m,t)= K(m,m')f(m',t)dm', (24) 

J mo 

C 2 (m,t)= m'K(m,m')f(m',t) dm', (25) 

J mo 

where the form of the kernel K(m,m') is the same that in § 2.2.1. We then fit Co and 
C 2 with a finite series in fractional powers of m in the same manner as given in Eq. (15). 
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Substituting these functions in place of one of the integrals (say over m'), we may integrate 
over m to get 



dM' Q 1 r h(i) 



dt 2M bv , IJlllu 

dM' 2 _ 1 f mL{t) 
~df~ M 2 (0) 



— — I f(m, t)C (m, t)dm = -~22 a^n M ' Pl 

rm L {t) 

j mf(m, t)C 2 (m, t)dm = ^2 hv Pl +iM' p 



Pi+1' 



r/l5(Pi-l)(P»-2) rn / f/l-P»(p»-2) f 7\/f'l ^Pi(Pi — 1) 



M' pi = [m^ m ' [M[]- p ^- Z} [M' 2 Y n ^- l} , (26) 

where we have made use of equations (8) and (9) to express the solution in terms of the 
integer moments k = 0,1,2. Here, [x Vi = M Pi (0)/M (0), u Pi+1 = M p . +1 (0)/M 2 (0), M' k = 
Mfc(t)/Mfc(0), and < Pi < 1. The Ck are fairly smooth functions over a large range of 
particle radii;, however, the accuracy of fitting a single series in fractional moments over 
a very broad range of particle sizes (i.e., over many orders of magnitude) may drop off 
significantly as the broadness of the range increases. This issue may be circumvented by 
employing a piecewise fit to the integrated kernels Ck- 

It is interesting to note that by this definition of Ck-, we have effectively accomplished 
what we set out to do in our discussion at the beginning of § 2.2, that is, defining the 
coagulation kernel in terms of finite series in powers of the mass m. The difference here is 
that we have done so through the first integral of the kernel, and not the kernel itself. This 
means that the mass distribution f(m, t) is expressed in the calulation of the Ck, but remains 
implicit in the definition of the ODEs (Eq. 26). Thus, the method is semi-implicit, because 
the RHS of the equations above can be expressed in terms of the moments (as defined in 
Eq. 2). Equations (26) are then integrated using the fourth order Runge-Kutta method, 
and may be compared with the results of § 2.2.1 and the direct integration of Eq. (1). 

This semi-implicit approach tracks the evolving kernel through the integration of Equa- 
tions (24) and (25) after every timestep, thus the computational time involved is similar to 
the explicit case. In order to update the kernel, one may solve for the new tjil after each At 
using the equation (q ^ 1) 



m L - m M q 



m 2 T - q - wt q (2 - q)M 1 



(27) 



The new normalization coefficient c can then be found from the definition of Mi = p. We 
then reintegrate Equations (24) and (25) under the powerlaw assumption, and then proceed 
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to fit the Ck{m,t) with a finite series in fractional powers of m. Although, in principle, any 
other two moments could be used to obtain thl, M q , which lies between M 2 and M 1; and 
because it roughly characterizes the evolution of the largest particle (see discussion at the 
end of § 2.2), seems the most consistent choice. The g-th moment is calculated using the 
Lagrange polynomial interpolation scheme (Eqs. 8 and 9). 

The advantage of the moments method lies in the ability to express the differential 
equations in terms of the moments of the distribution (i.e., their integrated properties). If a 
more explicitly mass-dependent approach is adopted (as is the case in § 2.2.1, and the semi- 
implicit approach described here), then the computational time significantly increases. One 
can improve the speed of computation by calculating the Ck periodically, or in the extreme 
case, only at t — which would make the approach truly implicit (e.g., § 2.1). The advantage 
of an implicit approach is that it becomes fully general (the form of / is only assumed at the 
onset), and also in the time it takes to solve (< 1 minute). The bulk of the time is spent in 
the integration of equations (24) and (25) which would occur only once. The disadvantage, 
of course, is that the particle velocity distribution is not updated as it changes with time 
(due to, e.g., changes in the bounds of the size distribution). We present examples of both 
extremes in § 3. 

If one wanted to implement a mass- or velocity-dependent sticking coefficient S, it can 
readily be included in the integration to obtain the Ck- The additional inclusion of source 
and sink terms due to erosion, fragmention, or gravitational growth in this semi-implicit 
formalism would require that we fit these terms in a similar manner to the Ck so that their 
subsequent integration over all m will yield sums over integer moments weighted by different 
sets of coefficients. Caution must be exercised in fitting, e.g., the gravitational growth term 
to ensure that the system of equations remains closed. Under such circumstances, a fully 
explicit approach such as that of § 2.2.1 may be preferred. Alternatively, these effects may 
be included in particle-histogram space in between coagulation interations. 

Finally, we should point out that allowing for other parameters (such as the index of 
the powerlaw q) to vary with time, does not pose a problem in either of the approaches we 
have presented. In both cases, one would simply need an additional moment, e.g. M 3 , to 
determine q(t). Similarly, a bifurcate d distribution in whic h the powerlaw exponent changes 



at a particular particle size (see, e.g.. lKenyon fc Luulll999l ) may also be studied. 



In § 3, we will compare the two approaches to the direct integration of the coagulation 
equation for cases in which there are only systematic velocities (vt = 0), as well as cases in 
which the velocity differences are induced by turbulent motions. 
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Numerical Results 



We carried out several calculations in order to demonstrate the accuracy of each alterna- 
tive method compared to the brute force integration of the collisional coagulation equation. 
For the purposes of comparison, unless otherwise noted, we chose the initial conditions to be 
a minimum mass solar nebula at 1 AU at a height of z — 10 3 km above the midplane, and 
a particle size distribution with a minimum initial radius of 1 cm and a maximum initial 
radius of 10 cm. The powerlaw exponent q = 11/6, which is assumed to be constant in 
these calculations, is representative of a fragmentation population. Standard integrations 
were carried out with a timestep of At = 10 years. 



3.1. Laminar Case (vt = 0) 

In Figure 2, we compare the case in which there are only systematic (pressure-gradient 
driven) velocities between particles (vt = 0). The solid curves represent the explicit assump- 
tion (invariant powerlaw slope q assumed in integration of Equations 22 and 23), while the 
dashed curves represent the implicit assumption. That is, in the integration of Eq. 26, the 
form of the mass distribution / is assumed only at t — 0. As before, the symbols represent 
the brute force calculation for two different grid resolutions (solid = 100 bins, open = 1000 
bins). Since with the explicit assumption we do not solve for the moments specifically, the 
values for the solid curves were obtained by substituting the time integrated values of c(t) 
and mi,{t) back into Eq. (17) for k = 0,2 only. This is because although we have used the 
moment equations to obtain these results, the explicit assumption of an invariant powerlaw 
size distribution means that only the equations for the moments Mq, Mi are needed (al- 
though see § 3.2). This is not true, however, for the implicit (and semi-implicit) assumption, 
where all the moments M , Mi, M 2 appear in the differential equations. 

We see that there is excellent agreement between both approaches and the numerical 
values obtained with the highest resolution case. In fact, the agreement between the explicit 
and implicit assumptions is also quite good. However, at this stage there has not been a 
great deal of growth (the largest particle size in the distribution has only grown to rj, = 11 
cm in size by the end of the simulation for the explicit case). Note that even though in 
the implicit case the Ck are calculated only once, the estimate for the largest mass tjil for 
the evolving distribution f(m,t) is still calculated using Eq. (27). At least in the case of 
minimal or slow growth, it appears that an implicit approach, or even periodic calculation 
of the Ck, may be sufficient. 

Numerical glitches in the low resolution brute force case arise from the interpolation 



-17- 



scheme for sampling the kernel. In particular, these glicthes are likely enhanced because 
the systematic relative velocities quickly approach zero for identical particle sizes (with the 
effect much more prominent for larger particles, hence not appearing so much in M ). The 
higher resolution case contains enough points to smooth out this effect. 

Note that (in all simulations) the direct integration of the coagulation equation gives a 
constant Mi as is to be expected given that we found dMi/dt = in the derivation of the 
moment equations in the absence of sources and sinks; this further validates the numerics of 
the brute force solution even for the complicated collisional kernel being utilized, and also 
provides a posteriori validation of our result that Mi = constant from, e.g., Eq. (5), which 
further validates derivation of equations (24) and (25) which was based on symmetry of the 
kernel. 



3.2. Turbulence Case (v T ^ 0) 

Next, we explored cases in which the systematic velocities were set to zero so that 
velocity differences between particles are due only to those induced by turbulence. Depending 
on the magnitude of the turbulence parameter a, these induced velocities can be either large 
or small relative to the systematic velocities. To demonstrate, we ran cases for three different 
values a = 10~ 6 , 10~ 5 , and 10~ 4 . In the absence of any mechanism to counter growth (e.g., 
fragmentation), larger a translates to faster rates of growth (due to larger relative velocities). 

In Figure 3 we plot the results for a = 10~ 4 , for the explicit (solid curves) and implicit 
(short dashed curves) approaches. The growth rate is more rapid than in the laminar case, 
with the second scaled moment ~ 70% larger (compared to Fig. 2) at the end of the run, 
meaning that the largest particle size achieved is r^ ~ 13.5 cm. We note that, as was the 
case for the M' 2 curves in Fig. 2, the explicit and implicit cases are very similar; however, 
this is not the case for the M' curves. The explicit approach overestimates the the value of 
Mq compared to the highest resolution brute force calculation, whereas the implicit approach 
significantly underestimates the zeroth moment. Both approaches understimate the value 
of M' 2 . When we compare integrations of M' 2 for smaller values of a as we do in Figure 4, 
we find that both approaches fit the coagulation calculation quite well. The growth rate for 
these two values of a are more in line with growth rate found when there are only systematic 
velocities, so the agreement should not be surprising. 

The long-dashed curves in Figures 3 and 4 are the implementation of the semi-implicit 
approach (§ 2.2.2) in which the Ck are updated at every time step. In this case, the semi- 
implicit approach provides a much a much better fit to the brute force calculation, indicating 
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that the kernel is evolving fast enough that an implicit approach cannot capture this effect. 
We consider the slight discrepencies between the semi-implicit approach and the brute force 
calculation to be as much a result of grid resolution as inaccuracy in using the Lagrange 
polynomial fits to the fractional moments. The explicit and fully implicit approach values 
of tl apparently understimate the largest particle size relative to the semi-implicit approach 
(which gives a value of Tl — 14 cm). 

Finally, we explored a variation of the explicit approach in which the condition Mi = p 
was strictly enforced (recall that Mi does not appear in Equations 18-19). This amounts to 
replacing Eq. (18) (derived for M ) with the corresponding equation for Mi. As a result, we 
cannot simultaneously fit both Mi and Mq. Alternatively, we can use only Eq. (19) for the 



second moment, by substituting c = (2 — q)p/(m 



2-q 
L 



m, 



2-<A 



(from the definition of Mi, Eq. 



17) in Eq. (19) so that Eq. (19) alone determines the growth of the largest particle, 
differential equation for mi then becomes 



The 



dvrii 
~dt 



(3-q)(2-q)pT 2 



(3 — q)(m 



2-q 
L 



m o 9 ) m l 9 ~ ( 2 "?)(% q 



m, 



3— q\ l—o 



(2* 



The results of the integration of Eq. (28) are shown as the dotted curves in Figure 3 and for 
the a = 10~ 4 case in Figure 4. It is clear that this modification provides a much better fit to 
M2, perhaps even better than the semi-implicit case above. However, using mi(t) calculated 
from Eq. (28), and solving for c(t) does a poor job fitting Mo (see Figure 3). Thus we 
conclude that although the modification of the explicit approach matches the brute force 
calculation of M 2 quite well, only the semi-implicit case is able to provide a simultaneous fit 
to both the zeroth and second moments (under the condition that Mi — p — constant). 



3.3. A Model Comparison 



Garaudl (120071 ) has developed a simplified analytical approach for dealing with the 



growth of particles in a turbulent regime, in which the particle size distribution is parame- 
terized by a powerlaw in particle mass with the same exponent we have used in this paper 
(q = 11 /6). Sim i lar to what we have presented in previous sections, the underlying assump- 
tion of iGaraudl (120071 ) is that collisions between particles occur frequently enough that a 
steady-state balance is reached in the form of the particle size dist r ibutio n, but with an 
upper size cutoff that varies with time. Thus, the model of IGaraudl (120071 ) is explicit and 
follows only two parameters: the growth of the largest particle mi(t), and the normalization 
factor c{m L ,t). 
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Given the similarity of the underlying assumptions for the Garaud model and the exam- 
ples we have presented, it is a useful exercise to compare the results of the two approaches 
directly The growth of the largest particle r L in the Garaud model (her Eq. 36), expressed 
in terms of the notation used in this paper, is given by 



dr L = pQH / «7St L 

dt ~ Ps y l + 64St|(2 + 5Sti a1 )- 2 ' [ ' 

In Eq. (29) Stx, is the Stokes number for the largest particle r L = (3m L /47rp s ) 1 / 3 , H = c g /Vt 
is the scale height of the gas, 7 is the adiabatic index of the gas, and in this expression it 
is assumed that the sticking coefficient S — 1, and t hat raj, » mo. The a bove equation is 



similar to the formula for grain growth proposed by IStepinski &: Valageasj (119971 ) to factors 
of order unity. Finally, we point out that the Garaud model for particle growth is restricted 
to the value q = 11/6 in order to preserve its completely analytical nature. As a means of 
a fair comparison, we chose to compare the Garaud model to our explicit case (Eq. 28) for 
reasons explained below. In the limit of mi >> mo and q = 11/6, Eq. (28) becomes 

drr » To , . 

— - ~ 0.05— — . 30 

at p s m L 

We present the results of our comparison in Figure 5 for the same initial conditions as 
described at the beginning of § 3.2 (upper curves). We find quite generally that our explicit 
calculation (Eqs. 28 and 30) leads to a faster growth rate than what is predicted from the 
Garaud model, initially. We note that the Garaud expression steepens quickly for r^ > 13 
cm, suggesting that the growth rates of the two approaches are more comparable at later 
times. However, the minor "kink" in the long-dashed curve is due to the shift from Epstein 
to Stokes flow. A much more subdued kink is visible in the explicit approach which uses 
the full expressions for the turbulent velocity, whereas in the derivation of Eq. (29), it is 
assumed that the stopping time in the Stokes regime is defined by some mean characteristic 
velocity which leads to a much more noticeable discontinuity. Regardless, the overall more 
subdued growth rate elicited by Eq. (29) (despite the steepening at later times due to a shift 
in flow regimes) is apparent from the lower set of curves in Fig. 5 where the initial conditions 
were chosen with r^(0) = 1 cm. For this case, growth occurs only in the Epstein regime, 
but still, the curves for the explicit approach and that calculated from Eq. (29) begin to 
diverge. Thus, the Garaud expression apparently underestimates the growth rate relative to 
our approach. 



We emphasize that the treatment of particle growth by I Gar audi ( 120071 ) requires several 



approximations in order to derive a purely analytical expressions for dr^/dt. Besides the 
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aforementioned restriction of q = 11/6, Garaud approximates the full expressions for the 
turbulent velocities that we use here by partitioning t> T into seperate cases dependent on 
the particles' stopping time relative to the turnover times of the smallest and largest scale 
eddies. Furthermore, some question may be raised as to the comparability of the moment 
equati ons used here to derive Eq. (28), versu s the particle g rowth equation used by iGaraud 
( 120071 . her Eq. 29). The equation used by I Gar audi (120071 ) is more akin to a "sweep-up" 
equation, with no sources or sinks, than to a formal coagulation equation. Because it bears 
some resemblance to the equation for M 2 (with some algebra, to factors of order unity), we 
concluded that our Eq. (28) is the appropriate analog. Despite the differences in growth 
rate, we find the agreement in the general trend of growth of the two approaches reassuring. 



4. Opacity Calculations 

Evolutionary models of protoplanetary nebu lae, giant planet atmospheres, etc. must 



somehow treat the escape of thermal radiation (jPollack et al.l 1 19961 ; iHubickyj et al.l 12005 



Durisen et al.ll2007l ). In most cases, particles provide the primary opacity for these mod- 
els. Observations of these and similar objects often rely on Spectral Energy Distributions 
(SEDs) which can be compared to a model once the model's internal temperature distribution 
is kn own; clear evidence is seen for grain growth in many cases (see review by iNatta et al. 
20071 ). Because of the nearly insurmountable computational burden involved with performing 
a fully self-consistent calculation of particle growth by coagulation along with an already dif- 
ficult fluid dynamical calculation, most modelers simply assume some invariant particle size 
distribution, such as th e MRN interstellar gr ain distribution, or make arbitrary assumptions 
about particle growth (IHubickyj et al.l 120051 ) . 



In the simplest regime (monodisperse particle radius r larger than a wavelength), the 
particle opacity can be written as the area per unit mass: 



4p s r 



cm 2 g x ; 



(31) 



thus growth in radius from O.l/i to 1 mm leads to a factor of 10 4 change in opacity. To 
the degree that this wavelength-independent regime holds, including particle size evolution 
by the moments method in one's evolutionary models would allow a very simple way to 
track particle growth and decreasing opacity. For instance, equation (31) above is easily 
generalized to the area per unit mass integrated over the size distribution: 
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/ ixr 2 f(m) dm _ / 9vr \ 1/3 M 2/3 



Jmf(m)dm \lQp 2 s ) M x ' ^ ^ 

As an application of the moments method in Figure 6, we have calculated the decrease in 
opacity (given by Eq. 32) with time using the semi-implicit approach. An initial particle size 
distribution with a lower bound of ro = 0.1 cm and q = 11/6 was used. Both the pressure 
gradient driven systematic velocities and the turbulence-induced velocities were used. In the 
absence of any mechanism to hinder particle growth, larger values of a lead to more steeply 
decreasing opacities with time. 

In a regime where the particle extinction efficiency Q(r,X) is wavelength-dependent 
(say, if the particles are comparable to or smaller than the wavelength), one simply inte- 
grates Q(r, A) over the powerlaw mass distributions resulting from the moments model. For 
example, 

Q7rr 2 Q(r,\)f(m)dm_ 1 /■•»* 2 



(33) 



«a = — °—fmi — ., , , = TT nr Q(r, X)c{m L )m q dm 

J mj{m) dm Mi 



mo 

1/3 r m L 



^(ST*^ 



These opacities K\ can be used to calculate Planck or Rosseland (wavelength-averaged) means 
for use in radiative transfer models. Recall that the powerlaw slope q can be freely adjusted 
within a small but plausible range to explore different growth regimes. 



5. Porosity 

Fractal growth of particles by low-velocity sticking of small solid monomers with radius 
r and mass m (and/or aggregates of such mon omers) causes them to have a density much 



less than the mate rial density of the monomers (IBeckwith et al.l |2000| ; iDominik et al.l 120071 ; 



Ormel et al.ll2008l ). These porous particles can be described as fractals with dimension D, 
such that the particle mass m increases proportionally to r D where r is some effective radius 
and D is the fractal dimension. Thus the particle's internal density is a function of particle 



size: 



/ \ ^ m <jm {r/r ) om \d-3 ^P° i i \D-3 /oa\ 

PV) = ~A 3 A 3 = A 3"( r / r °) = -A—^no) ■ ( 34 ) 

Anr 6 Airr 6 Aixr* Air 



-22- 



For a typical situation where D ~ 2 (JDominik et al.l 120071 ). p(r) oc r x and thus the product 



rp(r) is a constant across a wide range of particle sizes (until compaction sets in). These more 
complex but quite plausible particle density-size relationships complicate the expressions for 
particle stopping time and Stokes number (§ 2.2, Eq. 11; also, see appendix). 

Because the particle stopping times enter in through the collisional kernel, the fractal 
nature of particles can be accounted for using the method of moments in a straightforward 
manner while maintaining the criterion for closure of the system. We can verify this by 
noting that the mutual particle cross section 

a(m, m') = -^ fa 1 ' + m' 1/D ) 2 , (35) 



m 



2/D 



is proportional to m 2//D , whereas the stopping time (Eq. 11) in the Epstein regime (the regime 
that would apply to flu ffy fractal aggregates) where the drag force F D = (4/3)7ir 2 c g p g AV P g 



.e.g., 



Cuzzi et al.lll993l ) is given by 



t s = 3m = 3m ° m l - 2/D . (36) 

4:TTr 2 c g p g 4wr*CgPg 

In the limit of small Stokes number (St -C 1), both the systematic and turbulence-induced 
velocities are proportional to St, so that AV^ oc t s oc m 1_2//D , and the entire kernel is pro- 
portional to m. For larger St, the kernel has a shallower powerlaw dependence. Thus, in 
general, the dependence of K(m, m') on the mass is < m, preserving closure of the system, 
even for fractal particles. Even though the variation in the properties of the evolving dis- 
tribution due to porous particles are incorporated directly into the kernel, and are folded 
into the integration of the explicit approach, the effects of fractal aggregates can, nonethe- 
less, affect which moments characterize what properties in the semi-implicit (or implicit) 
approach. For example, the wavelength-independent opacity expressed in Eq. (32) takes the 
form k = (irr 2 /m 2 /D )M 2/D /M 1 . 

It should be noted that the value D = 2 represents a special case in that the particle- 
to-particle relative velocities in the Epstein regime do not depend on the mass of the fractal 
particle. Indeed, this would seem to indicate that fractal growth can proceed unabated with 
the corresponding stopping time of the fluffy aggregate remaining the same as that of a single 
monomer, which would have a significant effect on other particle properties. In particular, the 
wavelength- independent opacity for D = 2 is constant. However, impacts will ev entually lead 



to co mpaction or even fragmentation depending on the relative velocities (e.g.. lOrmel et al. 



20071 ). Both fractal grains and non- fractal particles in the same mass distribution can be 
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treated in the explicit approach without any modifications, while a piecewise fit to the 
integrated kernel Ck(m) (§ 2.2.2) can be used in order to account for the change in regimes 
in the semi-implicit approach. 



6. Conclusions 

We have demonstrated an approach to solving the collisional coagulation equation with 
an arbitrary collisional kernel which should be useful in cases when it is only necessary to 
keep track of general properties of the distribution. This approach involves solving a finite 
set of coupled differential equations in terms of the integer moments of the particle size 
distribution. The number of equations (and thus moments) needed depends on the number 
of properties being tracked. The advantage of the moments method approach is that it 
allows for considerable savings in computational time compared to direct integration of the 
coagulation equation, which requires keeping track of every particle size at every spatial 
location and timestep. 

In this paper we have specifically studied the growth of the largest particle under the 
assumption that the particle size distribution is a powerlaw; however, the technique can be 
extended to track other properties of the distribution that may change with time. There 
are many reasons to believe that a powerlaw size distribution is a natural end-state of par- 
ticle growt h, especially those with equal m ass per decade, because they have self-preserving 
properties (jCuzzi fc WeidenschiUindl2006l ). With the assumption of a powerlaw distribution, 



we have provided two different approaches to solving the moment equations, one explicit in 
which the powerlaw assumption is enforced rigorously at all times, and a semi-implicit ap- 
proach in which the kernel is integrated over one of the mass variables as much as once every 
time step. The latter approach can be made fully implicit by only assuming the form of the 
mass distribution at t — 0. These approaches are significantly faster than solutions of the 
coagulation equation because, in particular, the convolution integral (first term on RHS of 
Eq. 1) has been eliminated. In realistic evolutionary models, intermediate steps performed 
in particle "histogram" or "size-distribution" space may be interleaved with moments-based 
coagulation steps in order to account for, e.g., advective/transport terms. 

We have compared these alternate approaches to the brute force integration of the full 
coagulation equation for cases in which there are only systematic velocities, and cases in 
which the differences in velocity between particles is induced by turbulence. If the growth 
rate is gradual, the explicit and implicit approaches match the brute force calculation well 
(Fig. 2), whereas faster growth rates are more difficult to model (Fig. 3). We find that we 
are able to use the semi-implicit approach in which the Ck are updated at every timestep, 
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and a modification to the explicit approach, in which we solve the equation for M2 only with 
the assumption of p = Mi strictly enforced, in order to compensate for the faster growth 
rates. The modification to the explicit approach is useful if one is not particularly interested 
in following the evolution of the number density of particles, or other properties which may 
be decribed by (fractional) moments < 1 (e.g., see § 4). Our results also suggest that a fully 
implicit approach is probably most useful under circumstances in which the kernel depends 
on the mass in a straightforward manner (e.g, the Saffman- Turner kernel, § 2.1). 

We have compared the approach es developed in this paper to an alternative model for 



particle growth in a turbulent nebula (lGaraudll2007l ). We have found that there is fairly good 



agreement in general, but that the curves diverge as time proceeds. This does not appear to 
be particle-size dependent, or due to a shift in the flow regime. The Garaud expression (Eq. 
29) underestimates the growth rate of the largest particle size relative to our a pproach by only 



~ 20 — 30% in our comparison. We note that the advantage of the method of [Garaudj (120071 ) 
is that it is purely analytical; however, preservation of her analytical approach requires, 
amongst other things, the powerlaw exponent be restricted to q — 11/6. Our approach has 
no such restriction on the choice of exponent q, nor for q to even be a constant. 

As a sample application, we show how the moments method can be used in one's evo- 
lutionary model to track particle growth and opacity. As a specific case, we calculated the 
change in wavelength-independent opacity with time for an initial particle size distribution 
with upper and lower bounds of 0.1 — 1 cm. Both systematic and turbulent velocities were 
included. In the absence of any mechanism to counter growth, the opacity decreases sharply 
for higher choices of the turbulent parameter a. Extension to cases in which the extinction 
efficiency is wavelength dependent is straightforward. Such opacities can be used to calculate 
Planck or Rosseland wavelength-averaged means for use in radiative transfer codes. 

Finally, we indicate how porous particles with fractal dimension D can be accounted for 
in the moments method. The particle-size density relationships that arise affect the particle 
stopping times and Stokes number, both of which appear only in the collisional kernel. Thus 
implementation is straightforward. We can treat mass distributions composed of only fractal 
grains, or both fractal and non-fractal particles. In either case, relatively little modification 
is needed in the explicit approach, whereas with the semi-implicit (and implicit) approach, 
a piecewise fit to the integrated kernel Ck using the method as described in § 2.2.2 can be 
used to account for the change in particle growth regime when both types of particles are 
included. 

The computational burden of directly solving the coagulation equation makes it quite 
prohibitive to explore large regions of parameter space, and thus it serves as the primary 
bottleneck in evolutionary growth models. The approach demonstrated herein is intended 
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to obtain robust, quantitative results for disk properties such as particle growth timescales 
and "typical" particle sizes that may be used in modeling efforts that are focused more on 
the larger problem of planetesimal formation. 

Although we have not included sticking or fragmentation in this paper, we believe that 
the moment equations (Eq. 5) will remain valid at least up to the fragmentation barrier. 
This size will depend on the assumed particle strengths and choice of nebula parameters. 
The treatment of g rowth up to the fragmentation size is consistent with recent work by 



Brauer et al.l (120081 ) for the case in which additional effects such as radial drift are not 
included. The model as presented in this paper, however, should be useful for the early 
stages of protoplanetary nebula particle growth relevant to spectral energy distributions and 
MRI suppression. 

In a forthcoming paper we will explore the effects of a variable sticking coefficient 
S(m,m'). Furthermore, we will explore the addition of source and sink terms such as grav- 
itational growth, erosion, sublimation, and condensation in addition to fragmentation. The 
ultimate goal will be to apply this methodology to a global model that studies the evolution 
of both the gas and solids in nebular and subnebular (giant planetary) environments. 

We wish to thank Sandy Davis, Fred Ciesla, and Olenka Hubickyj for internal reviews 
which improved the exposition of this paper, and an anonymous reviewer for his or her careful 
analysis of the manuscript. This work was supported by a grant from NASA's Origins of 
Planetary Systems Program. 



A. Generalization of the Systematic (Pressure Gradient Driven) Velocities 



In this appendix, we generalize the basic equations of iNakagawa et al.l (119861 . also, see 
Tanaka et al. 2005) for a two-component fluid by extending the particle component to 
incorporate a size distribution. We adopt cylindrical coordinates (R,(f>,Z), designating the 
corresponding gas velocity components as (u, v, w) and the particle velocity components 
as (Ui,Vi,Wi). The equations of motion of particles and gas for the radial and tangential 
velocities, generalized for a particle size distribution, are given by 

dU- 

' -A iPg {Ui-u)+2W h (Al) 



dt 



-^ = -A iPg (Vi -v)- -QUi, (A2) 
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! = -5> ft („- ,7,) + «.-!§* (A3) 

^- = -Y^A jPj {v-V j )-\^ (A4) 

i 

where the system is assumed to be axysimmetric. Here, A t = (p 9 £ S j) _1 with t si the stopping 
time of a particle of radius r i: p g and p g are the gas mass density and pressure, and pj is 
the material density of particles of radius Tj. We have not included the vertical component 
equations; here, generalization to a particle size distribution is straightforward since it re- 
mains true in general that |iu/Wi-| << 1; that is, the v ertical gas component of the velocity 



w is negligibly small compared with Wi (see appendix. iNakagawa et al.lll986[ ). 



If the dust stopping time is short compared to the orbital period, we can seek steady- 
state solutions for the velocity components. Setting d/dt = 0, this system can be solved 
exactly in the Epstein regime. By solving Eq. (Al) for Vi and inserting into Eq. (A2), we 
obtain 

A 2 p 2 g u + 2ClAip g v u + 2^St, 
1 ~ A\pl + W ~ 1 + St 4 2 ' (A5J 

T/ _ -iyt/2)A lPg u + A 2 p 2 g v y - (l/2)nSt, 
1 " A}p 2 g + Q 2 " 1 + St, 2 ' (A6) 

where Vi was obtained by insertion of C/j back into Eq. (Al). In Eqs. (A5)-(A6) we have 
expressed the last equality in terms of the Stokes number Stj = t si Q, with Q being the local 
Kepler frequency. These expressions for the individual particle velocity components can be 
inserted into Eqs. (A3) and (A4) to yield expressions for radial and tangential gas velocity 
components. With little difficulty, one finds 

u = 2tjv k \ —, (A7) 



sl + (l + s ) 2 -' 

l + so 
's 2 + (l + s r 



V = -T)V K 2 rx, [At 



where 



^Eife^Ef^ («) 



A 2 p 2 + Q 2 ^ p g 1 + St 



j 



-27- 



^Eife^E^b- (mo) 



j- -jrg ■ j- rs- . "«j 



and 77 = — {l/2p g Vt R){dp g /dR). These expressions agree with those obtained by lTanaka et al 



( 120051 ) with the exception that their expressions for so and s\ are expressed in integral form. 



From these expressions, one can easily obtain expressions for the relative velocities with 
respect to the gas Ui — u and Vi — v , or the relative velocities between particles of different 
sizes Ui — Uj and V, — Vj. These expressions are also applicable for the Stokes regime in 
which the stopping times are a function of the relative particle-gas velocities (cf. Eq. 11). 
In this case, iterations must be performed in order to obtain the correct relative velocities. 
Convergence is generally achieved in a small number of iterations. 
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Fig. 1. — Comparison between the scaled integer moments M' k = Mk/Mk(0) for k — 0, 1,2 
(curves) obtained by the method of moments calculation with a simple Saffman- Turner tur- 
bulent coagulation kernel (Eq. 7) and a brute force integration (symbols) of the coagulation 
equation (Eq. 1). The integer moments for the latter are calculated a priori using the 
distribution function f(m,t). There is good agreement, especially for the higher resolution 



case. 



-31- 



0) 
-t-> 

c 

(D 

E 
o 



"D 

_(D 

O 

e> 



1.2 - 



0.96 



0.92 



0.88 - 



~\ 1 1 1 — r 



v T = 



~i 1 r 



explicit 

implicit 



• 1 2 bins 
o 1 3 bins 




10 



100 



time (years) 



Fig. 2. — Comparison of the scaled integer moments k = 0, 1, 2 obtained from the moments 
method (curves), and those obtained from the integration of the coagulation equation (Eq. 
1, symbols) for the case of a realistic collision kernel with vt = 0. Both the explicit (solid 
curves), and implicit (short-dashed curves) approaches match the coagulation calculation 
fairly well, especially the higher resolution case. 
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Fig. 3. — Comparison of the scaled integer moments M& for k — 0, 1, 2 obtained from the 
moments method, with results obtained from the brute force integration of the coagulation 
equation (Eq. 1) for the case of a realistic collision kernel with v t ^ (and no systematic 
velocities). In this case the turbulence parameter a = 10~ 4 . Both the explicit (solid curves) 
and implicit (short-dashed curves) approaches have some difficulty matching the coagulation 
calculation for both M and M 2 (note that the solid and short-dashed curves for M 2 lie on 
top of each other). However, the semi- implicit approach (long-dashed curves) provides a 
better fit. The modified explicit approach (dotted curves) provides a better fit for M 2 while 
giving a worse fit for M (see § 3.2). 
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Fig. 4. — Comparison of the second moment M' 2 (t) obtained from the moments method with 
that obtained from integration of the coagulation equation (Eq. 1) for the case of a realistic 
collision kernel with v t ^ (and no systematic velocities) for three different values of the 
turbulence parameter a. Both the explicit (solid curves) and implicit (short-dashed curves) 
approaches match the lower a values fairly well, but as indicated in fig. 3, have difficulty 
matching the case of a — 10 -4 (note that the solid and short-dashed curves lie on top of 
each other). However, not surprisingly, the semi-implicit approach (long-dashed curves) and 
a modified explicit approach (for M2 only, dotted curve) provides a much better fit (see § 
3.2). 
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Fig. 5. — Comparison of parti cle growth as a function of time using the explicit approach 
(Eq. 28, solid curve), and the iGaraudl (120071 ) analytical particle growth expression (long- 
dashed curve). The upper and lower sets of curves differ in the initial size of the largest 
particle. Both sets of curves begin to diverge immediately. The somewhat subdued kink in 
the upper curves at r L ~ 13.5 cm is due to a shift from Epstein to Stokes flow. However, 
overall, the agreement is not bad (~ 20 — 30% in r£). 
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Fig. 6. — Plot of the normalized (to initial value) wavelength-independent opacities for 
different choices of the turbulence parameter a. Opacities decrease sharply (in the absence 
of any mechanism to hinder particle growth) for higher a. Calculations were performed using 
the semi-implicit approach with both turbulent and systematic relative velocities included. 



